Supplementary Analyses

This document provides the code for the analyses reported in:

Cosme et al. (Preprint) Testing the adolescent social reorientation model using hierarchical growth curve modeling with parcellated fMRI data

This script reproduces the analyses reported in supplementary material.

load packages

library(tidyverse)
library(knitr)
library(lme4)
library(lmerTest)

define color palettes and labels

pal_self_other = c("#FFA90A", "#247BA0")
pal_social_academic = c("#63647E", "#F25F5C")
pal_wave = c("#693668", "#A74482", "#F84AA7")
pal_label = c("#47A8BD", "#DBC057", "#FF3366")
pal_gender = c("#70c1b3","#247BA0")

parcel_labeller = labeller(label = c('social' = 'social parcels', 'other' = 'control parcels', 'self' = 'self parcels'),
                           domain = c('social' = 'social domain', 'academic' = 'academic domain'),
                           wave = c("t1" = "wave 1", "t2" = "wave 2", "t3" = "wave 3"))

dcbw = theme_classic() +
  theme(text = element_text(size = 14, family = "Futura Medium", color = "black"),
        panel.background = element_blank(),
        plot.background = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 14),
        legend.background = element_rect(fill = NA, color = NA),
        axis.line = element_line(color = "black"),
        axis.text = element_text(color = "black"),
        panel.grid.minor = element_blank())

load MRI data

• exclude participants with >20% volumes with motion artifacts or low quality FX data
• standardize betas within parcel (across participants and time)
• recode standardized betas +/- 3 SDs from the mean as NA (0.8% of all observations)

# define parcels
self_parcels = c(5, 17, 28, 47, 52, 66, 116, 147, 152, 156, 184, 198, 207, 225, 249, 292, 309, 354, 380)
social_parcels = c(18, 23, 29, 49, 54, 59, 62, 63, 67, 76, 111, 114, 139, 143, 146, 150, 178, 179, 189, 203, 212, 224, 229, 236, 238, 239, 245, 250, 259, 266, 271, 301, 305, 310, 322, 324, 328, 331, 333, 342, 343, 350, 374, 391)

# mri exclusions
mri_exclusions = c('s002_t1', 's004_t1', 's008_t1', 's011_t1', 's017_t1', 
                   's026_t1', 's033_t2', 's034_t1', 's041_t1', 's044_t1', 
                   's047_t1', 's051_t1', 's054_t1', 's057_t1', 's059_t1', 
                   's061_t1', 's063_t1', 's070_t2', 's074_t1', 's074_t2', 
                   's078_t1', 's084_t1', 's090_t2', 's090_t3', 's094_t1', 
                   's094_t2', 's096_t1') 

# mri exclusions
parcellations = read_csv('../data/fxParcellations.csv') %>%
  mutate(label = ifelse(parcellation %in% self_parcels, 'self',
                 ifelse(parcellation %in% social_parcels, 'social', 'other'))) %>%
  mutate(wave = paste0("t", as.numeric(c(`10` = 1, `13` = 2, `16` = 3)[as.character(age)]))) %>%
  select(-age) %>%
  unite(sub_wave, c(subjectID, wave), remove = FALSE) %>%
  group_by(parcellation) %>%
  mutate(inclusion = ifelse(sub_wave %in% mri_exclusions, "excluded from MRI", "completed MRI"),
         beta = ifelse(sub_wave %in% mri_exclusions, NA, beta),
         sd = ifelse(sub_wave %in% mri_exclusions, NA, sd),
         beta_std = scale(beta, center = FALSE, scale = TRUE),
         mean_beta_std = mean(beta_std, na.rm = TRUE)) %>%
    select(-sub_wave) %>%
  ungroup()

# exclude parameter estimates 3 SD from the mean
parcellations_ex = parcellations %>%
  mutate(beta_std = ifelse(beta_std > mean_beta_std + 3 | beta_std < mean_beta_std - 3, NA, beta_std))

load cleaned SPPC and task data

# demographics
race_ethnicity = readxl::read_xlsx("../data/Ethnicity_w1.xlsx") %>%
  extract(SID, "subjectID", "L([0-9]{3})-.*") %>%
  mutate(subjectID = sprintf("s%s", subjectID),
         hispanic_latinx = ifelse(Ethnicity == 1, "yes", "no"),
         race = recode(Race, `1` = "White", `2` = "Multiracial", `3` = "Pacific Islander",
                       `4` = "Native American", `5` = "Black or African American")) %>%
  select(subjectID, hispanic_latinx, race)

demo = read.csv("../data/SFIC_age.pds.gender.csv") %>%
  rename("subjectID" = SID,
         "wave" = wavenum,
         "gender" = Gender) %>%
  mutate(subjectID = sprintf("s%03d", subjectID),
         wave = paste0("t", wave),
         age_c = age - 13,
         age_c2 = age_c^2,
         pdss_c = pdss - 3,
         pdss_c2 = pdss_c^2) %>%
  full_join(., race_ethnicity)

# self-report
sppc = readRDS("../data/sppc.RDS")

# task
task = readRDS("../data/task.RDS") %>%
  ungroup() %>%
  mutate(percent_stdgm = (task_percent - mean(task_percent, na.rm = TRUE)) / 10, #center at grand mean, in units of 10%
         percent_std50 = (task_percent - 50) / 10,
         percent_std = (task_percent - 50) / 10) %>% #center at 50, in units of 10%
  group_by(target, domain) %>%
  mutate(percent_stdc = (task_percent - mean(task_percent, na.rm = TRUE)) / 10) #center at condition mean, in units of 10%

# merge data
merged = parcellations_ex %>%
  full_join(., task, by = c("subjectID", "wave", "domain", "target")) %>%
  full_join(., sppc, by = c("subjectID", "wave", "domain")) %>%
  full_join(., demo, by = c("subjectID", "wave")) %>%
  mutate(inclusion = ifelse(is.na(inclusion), "didn't complete MRI", inclusion)) %>%
  filter(!(subjectID == "s086" & wave == "t3")) #no MRI, task, or self-report data was collected

# subset data for modeling
neuro_model_data_pdss = merged %>%
  filter(!is.na(beta)) %>%
  select(subjectID, wave, pdss, pdss_c, pdss_c2, target, domain, parcellation, label, beta, beta_std)

neuro_model_data_ex_pdss = neuro_model_data_pdss  %>%
  na.omit()

comp_model_data = merged %>%
  select(subjectID, wave, pdss, pdss_c, pdss_c2, age, age_c, age_c2, domain, competence) %>%
  unique() %>%
  na.omit()

imp_model_data = merged %>%
  select(subjectID, wave, pdss, pdss_c, pdss_c2, age, age_c, age_c2, domain, importance) %>%
  unique() %>%
  na.omit()

status_model_data = merged %>%
  select(subjectID, wave, pdss, pdss_c, pdss_c2, age, age_c, age_c2, target, domain, task_percent, contains("percent")) %>%
  unique() %>%
  na.omit()

# dummy code target and domain
neuro_model_data_dummy_pdss = neuro_model_data_ex_pdss %>%
  mutate(target = ifelse(target == "self", .5, -.5),
         domain = ifelse(domain == "social", .5, -.5))

sample descriptives

full sample

merged %>% 
  select(subjectID, wave, age, gender) %>%
  unique() %>%
  group_by(wave) %>%
  summarize(n = n(),
            meanAge = mean(age, na.rm = TRUE),
            sdAge = sd(age, na.rm = TRUE)) %>%
  kable(format = 'pandoc', digits = 2)

wave	n	meanAge	sdAge
t1	90	10.07	0.31
t2	57	13.04	0.31
t3	44	16.46	0.58

merged %>% 
  select(subjectID, wave, age, gender) %>%
  unique() %>%
  group_by(wave, gender) %>%
  summarize(n = n()) %>%
  spread(gender, n) %>%
  kable(format = 'pandoc', digits = 2)

wave	F	M
t1	45	45
t2	30	27
t3	25	19

merged %>% 
  select(subjectID, hispanic_latinx, race) %>%
  unique() %>%
  group_by(hispanic_latinx, race) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  mutate(total = sum(n),
         percent = (n/total) * 100) %>%
  select(-n, -total) %>%
  spread(hispanic_latinx, percent) %>%
  kable(format = 'pandoc', digits = 1)

race	no	yes
Black or African American	3.3	NA
Multiracial	11.1	5.6
Native American	4.4	1.1
Pacific Islander	4.4	1.1
White	43.3	25.6

broken down by inclusion

age_table = merged %>% 
  select(subjectID, inclusion, wave, age, pdss, gender) %>%
  unique() %>%
  group_by(inclusion, wave) %>%
  summarize(N = n(),
            Mage = mean(age, na.rm = TRUE),
            SDage = sd(age, na.rm = TRUE),
            Mpdss = mean(pdss, na.rm = TRUE),
            SDpdss = sd(pdss, na.rm = TRUE))

gender_table = merged %>% 
  select(subjectID, inclusion, wave, age, pds, gender) %>%
  unique() %>%
  group_by(inclusion, wave, gender) %>%
  summarize(N = n()) %>%
  spread(gender, N) %>%
  rename("Female" = F,
         "Male" = M)

merged %>% 
  select(subjectID, inclusion, hispanic_latinx, race) %>%
  unique() %>%
  group_by(inclusion, hispanic_latinx, race) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  mutate(total = sum(n),
         percent = (n/total) * 100) %>%
  select(-n, -total) %>%
  spread(hispanic_latinx, percent) %>%
  kable(format = 'pandoc', digits = 1)

inclusion	race	no	yes
completed MRI	Black or African American	0.8	NA
completed MRI	Multiracial	7.6	2.5
completed MRI	Native American	3.4	0.8
completed MRI	Pacific Islander	2.5	0.8
completed MRI	White	27.7	16.0
didn’t complete MRI	Black or African American	2.5	NA
didn’t complete MRI	Multiracial	1.7	NA
didn’t complete MRI	Native American	NA	0.8
didn’t complete MRI	Pacific Islander	0.8	NA
didn’t complete MRI	White	8.4	3.4
excluded from MRI	Multiracial	0.8	2.5
excluded from MRI	Native American	0.8	NA
excluded from MRI	Pacific Islander	0.8	0.8
excluded from MRI	White	11.8	2.5

age_table %>%
  left_join(gender_table) %>%
  extract(wave, "Wave", "t([1-3]{1})") %>%
  arrange(Wave) %>%
  select(Wave, inclusion, N, Mage, SDage, Mpdss, SDpdss, Female, Male) %>%
  rename("MRI status" = inclusion) %>%
  kable(format = 'pandoc', digits = 2)

Wave	MRI status	N	Mage	SDage	Mpdss	SDpdss	Female	Male
1	completed MRI	57	10.08	0.32	2.16	0.74	33	24
1	didn’t complete MRI	12	10.00	0.25	2.00	0.77	4	8
1	excluded from MRI	21	10.07	0.33	1.83	0.58	8	13
2	completed MRI	44	13.06	0.33	3.57	1.05	25	19
2	didn’t complete MRI	8	13.07	0.20	3.69	0.65	4	4
2	excluded from MRI	5	12.79	0.28	2.62	0.85	1	4
3	completed MRI	34	16.33	0.46	4.75	0.39	19	15
3	didn’t complete MRI	9	17.06	0.60	4.78	0.67	6	3
3	excluded from MRI	1	15.72	NA	3.50	NA	NA	1

age and puberty correlations

cor_fun = function(data) pmap(var.names, ~ cor.test(data[[.x]], data[[.y]])) %>% 
  map_df(broom::tidy) %>% 
  cbind(var.names, .)

age_pds = merged %>%
  select(subjectID, wave, gender, age, pdss) %>%
  gather(dev, value, age, pdss) %>%
  unite(wave_dev, dev, wave, sep = " ", remove = TRUE) %>%
  rownames_to_column() %>%
  spread(wave_dev, value) %>%
  group_by(subjectID) %>%
  fill(everything(), .direction = "up") %>%
  fill(everything(), .direction = "down") %>%
  select(-rowname) %>%
  unique() %>%
  ungroup()

var.names = tidystringdist::tidy_comb_all(names(select(age_pds, -c(subjectID, gender)))) %>%
  filter(!(grepl("age", V1) & grepl("age", V2)))

age_pds %>%
  nest(-c(gender)) %>%
  mutate(
    test = map(data, cor_fun)
  ) %>% 
  unnest(test, .drop = TRUE) %>%
  mutate(r = sprintf("%.02f", estimate), #sprintf("%.02f [%.02f, %.02f]", estimate, conf.low, conf.high)
         r = ifelse(p.value < .05 & p.value >= .01, paste0(r,"*"),
             ifelse(p.value < .01 & p.value >= .001, paste0(r,"**"),
             ifelse(p.value < .001, paste0(r,"***"), r))),
         r = gsub("0\\.", "\\.", r)) %>%
  select(gender, V1, V2, r) %>%
  rename("Sex" = gender) %>%
  unique() %>%
  mutate(V1 = gsub("age", "Age", V1),
         V2 = gsub("age", "Age", V2),
         V1 = gsub("pdss", "PDS-S", V1),
         V2 = gsub("pdss", "PDS-S", V2),
         V1 = gsub("t", "", V1),
         V2 = gsub("t", "", V2)) %>%
  spread(V2, r) %>%
  arrange(Sex) %>%
  filter(Sex == "F") %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "–-", .))) %>%
  kable(format = "pandoc")

Sex	V1	PDS-S 1	PDS-S 2	PDS-S 3
F	Age 1	.50***	.20	.01
F	Age 2	-.12	.40*	.30
F	Age 3	-.11	.02	.24
F	PDS-S 1	–-	.42*	-.08
F	PDS-S 2	–-	–-	.25

age_pds %>%
  nest(-c(gender)) %>%
  mutate(
    test = map(data, cor_fun)
  ) %>% 
  unnest(test, .drop = TRUE) %>%
  mutate(r = sprintf("%.02f", estimate), #sprintf("%.02f [%.02f, %.02f]", estimate, conf.low, conf.high)
         r = ifelse(p.value < .05 & p.value >= .01, paste0(r,"*"),
             ifelse(p.value < .01 & p.value >= .001, paste0(r,"**"),
             ifelse(p.value < .001, paste0(r,"***"), r))),
         r = gsub("0\\.", "\\.", r)) %>%
  select(gender, V1, V2, r) %>%
  rename("Sex" = gender) %>%
  unique() %>%
  mutate(V1 = gsub("age", "Age", V1),
         V2 = gsub("age", "Age", V2),
         V1 = gsub("pdss", "PDS-S", V1),
         V2 = gsub("pdss", "PDS-S", V2),
         V1 = gsub("t", "", V1),
         V2 = gsub("t", "", V2)) %>%
  spread(V1, r) %>%
  arrange(Sex) %>%
  filter(Sex == "M") %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "–-", .))) %>%
  kable(format = "pandoc")

Sex	V2	Age 1	Age 2	Age 3	PDS-S 1	PDS-S 2
M	PDS-S 1	-.32*	-.37	-.26	–-	–-
M	PDS-S 2	.22	.41*	.50*	.16	–-
M	PDS-S 3	-.24	.37	.41	.30	.62**

plot

full sample

merged %>% 
  select(subjectID, wave, gender, age, inclusion) %>%
  group_by(subjectID) %>%
  unique() %>%
  mutate(n.waves = n()) %>%
  ungroup() %>%
  arrange(n.waves, gender, age, inclusion) %>%
  mutate(order = row_number()) %>%
  ggplot(aes(age, reorder(subjectID, -order), color = gender)) +
    geom_line(size = .3) + 
    geom_point(aes(shape = inclusion), size = 2, fill = "white") +
    scale_color_manual(values = pal_gender, name = "", labels = c("female", "male")) +
    scale_shape_manual(name = "", values = c(19, 17, 21)) +
    labs(x = "\npubertal stage (PSD score)", y = "participant\n") +
    dcbw +
    theme(axis.line = element_line(colour = "black"),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = c(.75,.92),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm"))

MRI only

merged %>% 
  filter(!inclusion == "didn't complete MRI") %>%
  mutate(inclusion = ifelse(grepl("excluded", inclusion), "excluded", "included")) %>%
  select(subjectID, wave, gender, age, inclusion) %>%
  group_by(subjectID) %>%
  unique() %>%
  mutate(n.waves = n()) %>%
  ungroup() %>%
  arrange(n.waves, gender, age, inclusion) %>%
  mutate(order = row_number()) %>%
  ggplot(aes(age, reorder(subjectID, -order), color = gender)) +
    geom_line(size = .3) + 
    geom_point(aes(shape = inclusion), size = 2, fill = "white") +
    scale_color_manual(values = pal_gender, name = "", labels = c("female", "male")) +
    scale_shape_manual(name = "", values = c(21, 19)) +
    labs(x = "\npubertal stage (PSD score)", y = "participant\n") +
    dcbw +
    theme(axis.line = element_line(colour = "black"),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = c(.75,.92),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm"))

neural puberty models

visualize raw

Visualize the developmental trajectories using the raw data

hypothesis 0

domain_parc_plot_raw = neuro_model_data_pdss %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, y = beta_std_avg, color = domain)) + 
  geom_smooth(aes(group = interaction(parcellation, domain), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = domain), method = "lm", formula = y ~ poly(x, 2), size = 1.5) +
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_y_continuous(breaks = seq(-2, 2, 1)) + 
  coord_cartesian(ylim = c(-2, 2)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  theme_minimal(base_size = 12) +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

target_parc_plot_raw = neuro_model_data_pdss %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, y = beta_std_avg, color = target)) + 
  geom_smooth(aes(group = interaction(parcellation, target), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = target), method = "lm", formula = y ~ poly(x, 2), size = 1.5) +
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = pal_self_other) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_y_continuous(breaks = seq(-2, 2, 1)) + 
  coord_cartesian(ylim = c(-2, 2)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  theme_minimal(base_size = 12) +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h0_raw = cowplot::plot_grid(domain_parc_plot_raw, target_parc_plot_raw,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 1

domain_plot_raw = neuro_model_data_pdss %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, y = beta_std_avg, group = domain, color = domain)) +
  geom_smooth(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2), alpha = .2) + 
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.6, .65)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

soc_acad_plot_raw = neuro_model_data_pdss %>%
  group_by(subjectID, pdss, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, domain, beta_std_avg) %>%
  unique() %>%
  spread(domain, beta_std_avg) %>%
  mutate(avg_diff = social - academic) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = avg_diff)) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, color = pal_social_academic[2], fill = pal_social_academic[2], size = 1.5) + 
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.6, .65)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h1_raw = cowplot::plot_grid(domain_plot_raw, soc_acad_plot_raw,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 2

target_plot_raw = neuro_model_data_pdss %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, y = beta_std_avg, group = target, color = target)) +
  geom_smooth(aes(fill = target), method = 'lm', formula = y ~ poly(x,2), alpha = .2) + 
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = pal_self_other) +
  scale_y_continuous(breaks = round(seq(-.4, .4, .2), 1)) +
  coord_cartesian(ylim = c(-.5, .5)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

self_other_plot_raw = neuro_model_data_pdss %>%
  group_by(subjectID, pdss, label, target, parcellation) %>%
  mutate(beta_std_avf = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, target, beta_std_avf) %>%
  unique() %>%
  spread(target, beta_std_avf) %>%
  mutate(avg_diff = self - other) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = avg_diff)) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, color = pal_self_other[2], fill = pal_self_other[2], size = 1.5) + 
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.4, .4, .2), 1)) +
  coord_cartesian(ylim = c(-.5, .5)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h2_raw = cowplot::plot_grid(target_plot_raw, self_other_plot_raw,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 3

int_plot_raw = neuro_model_data_pdss %>%
  distinct(parcellation, target, domain, pdss, beta_std, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = beta_std, group = interaction(target, domain), color = domain, linetype = target)) +
  geom_smooth(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2), alpha = .2) + 
  scale_y_continuous(breaks = c(-.4, .2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .45)) +
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_linetype_manual(name = "", values = c("dotted", "solid")) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\npubertal stage (PSD score)', y = 'mean standardized parameter estimate\n', color = '', linetype = '', fill = '') +
  guides(linetype = guide_legend(override.aes = list(color = "black"))) +
  dcbw +
  theme(legend.position = c(.75, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(.5, "cm"),
        legend.direction = "vertical",
        legend.box = "horizontal")

soc_acad_int_plot_raw = neuro_model_data_pdss %>%
  group_by(subjectID, pdss, label, target, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, target, domain, beta_std_avg) %>%
  unique() %>%
  spread(domain, beta_std_avg) %>%
  mutate(avg_diff = social - academic) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = avg_diff, color = target)) +
  geom_smooth(aes(group = interaction(parcellation, target), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = target), method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, size = 1.5) + 
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = pal_self_other) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.65, .65)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

self_other_int_plot_raw = neuro_model_data_pdss %>%
  group_by(subjectID, pdss, label, domain, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, domain, target, beta_std_avg) %>%
  unique() %>%
  spread(target, beta_std_avg) %>%
  mutate(avg_diff = self - other) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = avg_diff, color = domain)) +
  geom_smooth(aes(group = interaction(parcellation, domain), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, size = 1.5) + 
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.65, .65)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h3_raw = cowplot::plot_grid(int_plot_raw, soc_acad_int_plot_raw, self_other_int_plot_raw,
                       labels = c('A', 'B', 'C'), ncol = 3))

run models

Estimate full models

domain x puberty

model_domain_equation = formula(beta_std ~ 1 + pdss_c*domain*label +
                                  pdss_c2*domain*label +
                                  (1 + pdss_c*domain || subjectID) +
                                  (1 + pdss_c*domain + pdss_c2*domain || parcellation))

model_domain_formula = lFormula(model_domain_equation, data = neuro_model_data_dummy_pdss)

model_domain_numFx = length(dimnames(model_domain_formula$X)[[2]])

model_domain_numRx = sum(as.numeric(lapply(model_domain_formula$reTrms$cnms, function(x) {
  l <- length(x)
  (l*(l - 1)) / 2 + l
})))

model_domain_maxfun = 10*(model_domain_numFx + model_domain_numRx + 1)^2

if (file.exists("../data/model_domain_pdss.RDS")) {
  model_domain = readRDS("../data/model_domain_pdss.RDS")
} else {
  model_domain = lmer(model_domain_equation, data = neuro_model_data_dummy_pdss, REML = F, #Use ML since we want to compare random effects
                      verbose = 2,
                      control = lmerControl(optCtrl = list(maxfun = model_domain_maxfun), optimizer = "bobyqa", calc.derivs = FALSE))
  saveRDS(model_domain, "../data/model_domain_pdss.RDS")
}

domain x target x puberty

model_target_equation = formula(beta_std ~ 1 + pdss_c*target*domain*label +
                                  pdss_c2*target*domain*label +
                                  (1 + pdss_c*target*domain || subjectID) +
                                  (1 + pdss_c*target*domain + pdss_c2*target*domain || parcellation))

model_target_formula = lFormula(model_target_equation, data = neuro_model_data_dummy_pdss)

model_target_numFx = length(dimnames(model_target_formula$X)[[2]])

model_target_numRx = sum(as.numeric(lapply(model_target_formula$reTrms$cnms, function(x) {
  l <- length(x)
  (l*(l - 1)) / 2 + l
})))

model_target_maxfun = 10*(model_target_numFx + model_target_numRx + 1)^2

if (file.exists("../data/model_target_pdss.RDS")) {
  model_target = readRDS("../data/model_target_pdss.RDS")
} else {
  model_target = lmer(model_target_equation, data = neuro_model_data_dummy_pdss, REML = F, #Use ML since we want to compare random effects
                      verbose = 2,
                      control = lmerControl(optCtrl = list(maxfun = model_target_maxfun), optimizer = "bobyqa", calc.derivs = FALSE))
  saveRDS(model_target, "../data/model_target_pdss.RDS")
}

compare

anova(model_domain, model_target) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = Df,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")

model	npar	AIC	X2	X2 df	p
model_domain	29	438603.3	–	–	–
model_target	57	438378.3	280.92	28	< .001

summarize best fitting model

model_target %>%
  broom.mixed::tidy(effects = c("ran_pars", "fixed"), conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  select(-group) %>%
  rename("b" = estimate,
         "SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept (puberty 3, label (control))", term),
         term = gsub("target", "Target", term),
         term = gsub("domain", "Domain", term),
         term = gsub("labelself", "Label (self)", term),
         term = gsub("labelsocial", "Label (social)", term),
         term = gsub("pdss_c", "Puberty", term),
         term = gsub(":", " x ", term),
         term = gsub("sd__", "", term),
         term = gsub("Observation", "observation", term),
         effect = gsub("ran_pars", "random", effect),
         `b [95% CI]` = ifelse(effect == "fixed", sprintf("%.3f [%.3f, %.3f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 3) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  select(term, `b [95% CI]`, SE, t, df, p) %>%
  kable(format = "pandoc")

term	b [95% CI]	SE	t	df	p
Intercept (puberty 3, label (control))	-0.008 [-0.060, 0.045]	0.027	-0.280	398.452	.779
Puberty	0.006 [-0.005, 0.017]	0.005	1.122	114.748	.264
Target	-0.012 [-0.030, 0.006]	0.009	-1.343	163.818	.181
Domain	0.004 [-0.021, 0.030]	0.013	0.342	126.671	.733
Label (self)	-0.010 [-0.214, 0.194]	0.104	-0.097	356.861	.923
Label (social)	0.209 [0.069, 0.348]	0.071	2.936	356.762	.004
Puberty2	-0.000 [-0.006, 0.006]	0.003	-0.013	699.220	.990
Puberty x Target	0.004 [-0.006, 0.014]	0.005	0.833	87.273	.407
Puberty x Domain	0.004 [-0.013, 0.021]	0.009	0.478	65.176	.634
Target x Domain	-0.009 [-0.046, 0.027]	0.019	-0.506	146.016	.614
Puberty x Label (self)	0.004 [-0.021, 0.029]	0.013	0.308	378.167	.759
Puberty x Label (social)	0.015 [-0.002, 0.033]	0.009	1.705	377.575	.089
Target x Label (self)	0.214 [0.157, 0.271]	0.029	7.340	2050.846	< .001
Target x Label (social)	-0.019 [-0.058, 0.020]	0.020	-0.949	2021.793	.343
Domain x Label (self)	0.126 [0.058, 0.194]	0.035	3.612	692.867	< .001
Domain x Label (social)	0.144 [0.097, 0.191]	0.024	6.060	685.991	< .001
Target x Puberty2	-0.002 [-0.010, 0.006]	0.004	-0.508	360.135	.612
Domain x Puberty2	-0.012 [-0.022, -0.002]	0.005	-2.255	558.945	.025
Label (self) x Puberty2	-0.012 [-0.029, 0.004]	0.008	-1.480	383.964	.140
Label (social) x Puberty2	-0.001 [-0.012, 0.010]	0.006	-0.174	380.939	.862
Puberty x Target x Domain	-0.019 [-0.041, 0.004]	0.012	-1.610	76.241	.112
Puberty x Target x Label (self)	0.027 [-0.002, 0.055]	0.015	1.816	176512.793	.069
Puberty x Target x Label (social)	0.008 [-0.011, 0.028]	0.010	0.839	176505.823	.402
Puberty x Domain x Label (self)	0.021 [-0.007, 0.050]	0.015	1.470	176512.731	.142
Puberty x Domain x Label (social)	0.007 [-0.013, 0.027]	0.010	0.699	176504.705	.484
Target x Domain x Label (self)	-0.128 [-0.239, -0.018]	0.056	-2.280	176522.246	.023
Target x Domain x Label (social)	0.055 [-0.020, 0.130]	0.038	1.432	176510.597	.152
Target x Domain x Puberty2	0.005 [-0.012, 0.022]	0.009	0.600	381.300	.549
Target x Label (self) x Puberty2	-0.022 [-0.047, 0.004]	0.013	-1.647	176520.765	.100
Target x Label (social) x Puberty2	0.002 [-0.016, 0.019]	0.009	0.186	176510.160	.852
Domain x Label (self) x Puberty2	-0.003 [-0.029, 0.023]	0.013	-0.204	1002.405	.838
Domain x Label (social) x Puberty2	0.005 [-0.013, 0.022]	0.009	0.512	989.901	.608
Puberty x Target x Domain x Label (self)	-0.041 [-0.099, 0.016]	0.029	-1.414	176512.990	.157
Puberty x Target x Domain x Label (social)	-0.025 [-0.064, 0.014]	0.020	-1.263	176504.330	.207
Target x Domain x Label (self) x Puberty2	0.101 [0.050, 0.152]	0.026	3.849	176520.602	< .001
Target x Domain x Label (social) x Puberty2	0.020 [-0.015, 0.055]	0.018	1.123	176506.142	.261

# summarize random effects
print(VarCorr(model_target), comp = c("Variance","Std.Dev."), digits = 3)

##  Groups          Name                  Variance  Std.Dev.
##  parcellation    target:domain:pdss_c2 0.0000000 0.00000 
##  parcellation.1  pdss_c:target:domain  0.0000000 0.00000 
##  parcellation.2  domain:pdss_c2        0.0000992 0.00996 
##  parcellation.3  target:pdss_c2        0.0000000 0.00000 
##  parcellation.4  target:domain         0.0000000 0.00000 
##  parcellation.5  pdss_c:domain         0.0000000 0.00000 
##  parcellation.6  pdss_c:target         0.0000000 0.00000 
##  parcellation.7  pdss_c2               0.0004905 0.02215 
##  parcellation.8  domain                0.0075317 0.08679 
##  parcellation.9  target                0.0010029 0.03167 
##  parcellation.10 pdss_c                0.0020591 0.04538 
##  parcellation.11 (Intercept)           0.1899332 0.43581 
##  subjectID       pdss_c:target:domain  0.0029573 0.05438 
##  subjectID.1     target:domain         0.0054222 0.07364 
##  subjectID.2     pdss_c:domain         0.0027859 0.05278 
##  subjectID.3     pdss_c:target         0.0003843 0.01960 
##  subjectID.4     domain                0.0039597 0.06293 
##  subjectID.5     target                0.0010293 0.03208 
##  subjectID.6     pdss_c                0.0006930 0.02633 
##  subjectID.7     (Intercept)           0.0024120 0.04911 
##  Residual                              0.6669937 0.81670

visualize fitted models

Visualize the developmental trajectory using the fitted values from the domain x target x age model

reForm = as.formula("~(1 + pdss_c*target*domain + pdss_c2*target*domain || parcellation)")

neuro_plot_data_puberty = with(neuro_model_data_dummy_pdss, 
                    expand.grid(target = unique(target), 
                                domain = unique(domain),
                                parcellation = unique(parcellation),
                                pdss = unique(pdss),
                                stringsAsFactors = F)) %>%
  mutate(label = ifelse(parcellation %in% self_parcels, 'self',
                 ifelse(parcellation %in% social_parcels, 'social', 'other')),
         pdss_c = pdss - 3, 
         pdss_c2 = pdss_c^2, 
         subjectID = NA)

neuro_plot_data_puberty$expected = predict(model_target, newdata = neuro_plot_data_puberty, re.form = reForm)
neuro_plot_data_puberty$expected_mean = predict(model_target, newdata = neuro_plot_data_puberty, re.form = NA)

neuro_plot_data_puberty = neuro_plot_data_puberty %>%
  mutate(target = factor(target, levels = c(-.5, .5), labels = c("other", "self")),
         domain = factor(domain, levels = c(-.5, .5), labels = c("academic", "social")))

hypothesis 0

domain_parc_plot = neuro_plot_data_puberty %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, pdss, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, y = expected_avg, color = domain)) + 
  geom_smooth(aes(group = interaction(parcellation, domain), size = label), method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1.25, se = FALSE) +
  scale_color_manual(name = "", values = pal_social_academic) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n") + 
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

target_parc_plot = neuro_plot_data_puberty %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  group_by(subjectID, pdss, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  ggplot(aes(x = pdss, y = expected_avg, color = target)) + 
  geom_smooth(aes(group = interaction(parcellation, target), size = label), method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1.25, se = FALSE) +
  scale_color_manual(values = pal_self_other) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n", color = "") + 
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h0_fitted = cowplot::plot_grid(domain_parc_plot, target_parc_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 1

domain_plot = neuro_plot_data_puberty %>%
  group_by(subjectID, pdss, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected_mean, na.rm = TRUE)) %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_avg, color = domain)) +
  #geom_point(aes(group = target)) + 
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = seq(-.2, .45, .2)) +
  coord_cartesian(ylim = c(-.25, .4)) +
  scale_color_manual(values = pal_social_academic) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

soc_acad_plot = neuro_plot_data_puberty %>%
  group_by(subjectID, pdss, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, domain, expected_avg) %>%
  unique() %>%
  spread(domain, expected_avg) %>%
  mutate(expected_diff = social - academic) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_diff)) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, color = pal_social_academic[2], size = 1.5) + 
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = seq(-.2, .45, .2)) +
  coord_cartesian(ylim = c(-.25, .4)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = "none")

(h1_fitted = cowplot::plot_grid(domain_plot, soc_acad_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 2

target_plot = neuro_plot_data_puberty %>%
  group_by(subjectID, pdss, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected_mean, na.rm = TRUE)) %>%
  distinct(parcellation, target, target, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_avg, color = target)) +
  #geom_point(aes(group = target)) + 
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = seq(-.1, .3, .1)) +
  coord_cartesian(ylim = c(-.15, .3)) +
  scale_color_manual(values = pal_self_other) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

self_other_plot = neuro_plot_data_puberty %>%
  group_by(subjectID, pdss, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, target, expected_avg) %>%
  unique() %>%
  spread(target, expected_avg) %>%
  mutate(expected_diff = self - other) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_diff)) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, color = pal_self_other[2], size = 1.5) + 
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = seq(-.1, .3, .1)) +
  coord_cartesian(ylim = c(-.15, .3)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = "none")

(h2_fitted = cowplot::plot_grid(target_plot, self_other_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 3

int_plot = neuro_plot_data_puberty %>%
  distinct(parcellation, target, domain, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_mean, group = interaction(target, domain), color = domain, linetype = target)) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = 1, se = FALSE, size = 1.25) + 
  #geom_point() + 
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .45)) +
  scale_color_manual(values = pal_social_academic) +
  scale_linetype_manual(name = "", values = c("dotted", "solid")) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n",  color = '', linetype = '') +
  guides(linetype = guide_legend(override.aes = list(color = "black"))) +
  dcbw +
  theme(legend.position = c(.75, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(.5, "cm"),
        legend.direction = "vertical",
        legend.box = "horizontal")

soc_acad_int_plot = neuro_plot_data_puberty %>%
  group_by(subjectID, pdss, label, target, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, target, domain, expected_avg) %>%
  unique() %>%
  spread(domain, expected_avg) %>%
  mutate(expected_diff = social - academic) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_diff, color = target)) +
  geom_smooth(aes(group = interaction(parcellation, target), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, size = 1.5) + 
  scale_color_manual(values = pal_self_other) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .45)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))


self_other_int_plot = neuro_plot_data_puberty %>%
  group_by(subjectID, pdss, label, domain, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, pdss, label, domain, target, expected_avg) %>%
  unique() %>%
  spread(target, expected_avg) %>%
  mutate(expected_diff = self - other) %>%
  distinct(parcellation, pdss, label, .keep_all = T) %>%
  ggplot(aes(x = pdss, y = expected_diff, color = domain)) +
  geom_smooth(aes(group = interaction(parcellation, domain), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, size = 1.5) + 
  scale_color_manual(values = pal_social_academic) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .45)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\npubertal stage (PSD score)", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h3_fitted = cowplot::plot_grid(int_plot, soc_acad_int_plot, self_other_int_plot,
                       labels = c('A', 'B', 'C'), ncol = 3))

self/other evaluation task analyses

visualize raw

age

(status_plot_raw_age = status_model_data %>%
  filter(!is.na(task_percent)) %>%
  ggplot(aes(age, task_percent, color = domain, fill = domain, linetype = target)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_x_continuous(limits = c(min(status_model_data$age), 18), breaks = c(10, 12, 14, 16, 18)) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_fill_manual(name = "", values = pal_social_academic) +
    scale_linetype_manual(name = "", values = c("dotted", "solid")) +
    guides(linetype = guide_legend(override.aes = list(color = "black"))) +
    labs(x = "\npubertal stage (PSD score)", y = "social status or academic proficiency\n(% endorsed)\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

puberty

(status_plot_raw = status_model_data %>%
  filter(!is.na(task_percent)) %>%
  ggplot(aes(pdss, task_percent, color = domain, fill = domain, linetype = target)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_fill_manual(name = "", values = pal_social_academic) +
    scale_linetype_manual(name = "", values = c("dotted", "solid")) +
    guides(linetype = guide_legend(override.aes = list(color = "black"))) +
    labs(x = "\npubertal stage (PDS score)", y = "social status or academic proficiency\n(% endorsed)\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

run models

ICC

model.icc = anova(lm(task_percent ~ subjectID, data = status_model_data))
round(model.icc$`Sum Sq`[1] / sum(model.icc$`Sum Sq`), 2)

## [1] 0.18

age

`status model 1` = lmer(task_percent ~ 1 + target*domain + 
                        (1 + target * domain | subjectID), 
                        data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 2 age` = lmer(task_percent ~ 1 + target*domain + age_c + 
                        (1 + target * domain | subjectID), 
                        data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 3 age` = lmer(task_percent ~ 1 + target*domain + age_c + age_c2 + 
                         (1 + target * domain | subjectID), 
                         data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 4 age` = lmer(task_percent ~ 1 + target*domain*age_c + 
                            (1 + target * domain | subjectID), 
                           data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 5 age` = lmer(task_percent ~ 1 + target*domain*age_c + 
                            target*domain*age_c2 + 
                            (1 + target * domain | subjectID), 
                            data = status_model_data, control = lmerControl(optimizer = "bobyqa"))

puberty

`status model 2 puberty` = lmer(task_percent ~ 1 + target*domain + pdss_c + 
                        (1 + target * domain | subjectID), 
                        data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 3 puberty` = lmer(task_percent ~ 1 + target*domain + pdss_c + pdss_c2 + 
                         (1 + target * domain | subjectID), 
                         data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 4 puberty` = lmer(task_percent ~ 1 + target*domain*pdss_c + 
                            (1 + target * domain | subjectID), 
                           data = status_model_data, control = lmerControl(optimizer = "bobyqa"))
`status model 5 puberty` = lmer(task_percent ~ 1 + target*domain*pdss_c + target*domain*pdss_c2 + 
                            (1 + target * domain | subjectID), 
                            data = status_model_data, control = lmerControl(optimizer = "bobyqa"))

compare models

age

anova(`status model 1`, `status model 2 age`, `status model 3 age`, `status model 4 age`, `status model 5 age`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = Df,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")

model	npar	AIC	X2	X2 df	p
status model 1	15	5237.97	–	–	–
status model 2 age	16	5239.91	0.05	1	.815
status model 3 age	17	5241.69	0.22	1	.637
status model 4 age	19	5233.46	12.23	2	.002
status model 5 age	23	5231.07	10.39	4	.034

puberty

anova(`status model 1`, `status model 2 puberty`, `status model 3 puberty`, `status model 4 puberty`, `status model 5 puberty`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = Df,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")

model	npar	AIC	X2	X2 df	p
status model 1	15	5237.97	–	–	–
status model 2 puberty	16	5238.96	1	1	.317
status model 3 puberty	17	5238.77	2.2	1	.138
status model 4 puberty	19	5232.81	9.95	2	.007
status model 5 puberty	23	5237.02	3.8	4	.434

summarize best fitting models

age

`status model 5 age` %>%
  broom.mixed::tidy(effects = c("ran_pars", "fixed"), conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  select(-group) %>%
  rename("b" = estimate,
         "SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p)),
         term = gsub("\\(Intercept\\)", "Intercept (academic, other, age 13)", term),
         term = gsub("domainsocial", "Domain (social)", term),
         term = gsub("targetself", "Target (self)", term),
         term = gsub("age_c", "Age", term),
         term = gsub(":", " x ", term),
         term = gsub("sd__", "", term),
         term = gsub("Observation", "observation", term),
         effect = gsub("ran_pars", "random", effect),
         `b [95% CI]` = ifelse(effect == "fixed", sprintf("%.2f [%.2f, %.2f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 2) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  select(term, `b [95% CI]`, SE, t, df, p) %>%
  kable(format = "pandoc")

term	b [95% CI]	SE	t	df	p
Intercept (academic, other, age 13)	57.87 [51.28, 64.47]	3.36	17.21	174.40	< .001
Target (self)	4.82 [-4.34, 13.97]	4.67	1.03	182.28	.304
Domain (social)	2.29 [-6.42, 11.01]	4.45	0.52	188.62	.607
Age	-0.56 [-1.65, 0.53]	0.55	-1.01	365.43	.312
Age2	0.14 [-0.42, 0.69]	0.28	0.48	347.30	.631
Target (self) x Domain (social)	18.51 [6.41, 30.61]	6.17	3.00	194.63	.003
Target (self) x Age	-0.68 [-2.22, 0.86]	0.79	-0.87	366.83	.386
Domain (social) x Age	0.55 [-0.97, 2.07]	0.78	0.70	371.77	.481
Target (self) x Age2	0.43 [-0.35, 1.21]	0.40	1.07	348.69	.283
Domain (social) x Age2	-0.48 [-1.26, 0.29]	0.40	-1.22	352.91	.222
Target (self) x Domain (social) x Age	2.09 [-0.04, 4.22]	1.09	1.92	372.86	.056
Target (self) x Domain (social) x Age2	-0.64 [-1.73, 0.46]	0.56	-1.14	356.85	.255

print(VarCorr(`status model 5 age`), comp = c("Variance","Std.Dev."), digits = 5)

##  Groups    Name                    Variance Std.Dev. Corr                
##  subjectID (Intercept)              377.05  19.418                       
##            targetself               689.50  26.258   -0.815              
##            domainsocial             565.08  23.771   -0.866  0.798       
##            targetself:domainsocial 1040.78  32.261    0.772 -0.936 -0.905
##  Residual                           239.13  15.464

puberty

`status model 4 puberty` %>%
  broom.mixed::tidy(effects = c("ran_pars", "fixed"), conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  select(-group) %>%
  rename("b" = estimate,
         "SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p)),
         term = gsub("\\(Intercept\\)", "Intercept (academic, other, stage 3)", term),
         term = gsub("domainsocial", "Domain (social)", term),
         term = gsub("targetself", "Target (self)", term),
         term = gsub("pdss_c", "Puberty", term),
         term = gsub(":", " x ", term),
         term = gsub("sd__", "", term),
         term = gsub("Observation", "observation", term),
         effect = gsub("ran_pars", "random", effect),
         `b [95% CI]` = ifelse(effect == "fixed", sprintf("%.2f [%.2f, %.2f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 2) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  select(term, `b [95% CI]`, SE, t, df, p) %>%
  kable(format = "pandoc")

term	b [95% CI]	SE	t	df	p
Intercept (academic, other, stage 3)	59.23 [54.19, 64.28]	2.57	23.01	78.80	< .001
Target (self)	8.48 [1.51, 15.45]	3.56	2.38	80.81	.02
Domain (social)	-1.54 [-8.03, 4.94]	3.31	-0.47	80.29	.642
Puberty	-1.68 [-3.91, 0.54]	1.13	-1.48	380.92	.139
Target (self) x Domain (social)	12.33 [3.43, 21.24]	4.54	2.71	80.80	.008
Target (self) x Puberty	-1.14 [-4.28, 2.00]	1.60	-0.71	379.97	.478
Domain (social) x Puberty	1.37 [-1.72, 4.46]	1.57	0.87	382.90	.385
Target (self) x Domain (social) x Puberty	3.53 [-0.78, 7.83]	2.20	1.60	383.96	.109

print(VarCorr(`status model 4 puberty`), comp = c("Variance","Std.Dev."), digits = 5)

##  Groups    Name                    Variance Std.Dev. Corr                
##  subjectID (Intercept)              367.36  19.167                       
##            targetself               687.68  26.224   -0.825              
##            domainsocial             561.36  23.693   -0.869  0.796       
##            targetself:domainsocial 1030.63  32.103    0.797 -0.941 -0.908
##  Residual                           244.39  15.633

visualize fitted

age

model_to_plot = `status model 5 age`
REFormString = as.character(findbars(model_to_plot@call$formula)[[1]])
reForm = as.formula(paste0('~(', REFormString[[2]], REFormString[[1]], REFormString[[3]], ')'))
status_model_data$expected_mean = predict(model_to_plot, newdata = status_model_data, re.form = NA)

(status_plot_age = status_model_data %>%
  filter(!is.na(target)) %>%
  ggplot(aes(age, expected_mean, color = domain, linetype = target)) +
    geom_line(aes(age, task_percent, group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_point(aes(age, task_percent, group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2, se = FALSE, size = 1.25) +
    scale_x_continuous(limits = c(min(status_model_data$age), 18), breaks = c(10, 12, 14, 16, 18)) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_linetype_manual(name = "", values = c("dotted", "solid")) +
    coord_cartesian(ylim = c(0, 100)) +
    guides(linetype = guide_legend(override.aes = list(color = "black"))) +
    labs(x = "\npubertal stage (PSD score)", y = "social status or academic proficiency\n") +
    dcbw +
    theme(legend.position = "none"))

puberty

model_to_plot = `status model 4 puberty`
REFormString = as.character(findbars(model_to_plot@call$formula)[[1]])
reForm = as.formula(paste0('~(', REFormString[[2]], REFormString[[1]], REFormString[[3]], ')'))
status_model_data$expected_mean = predict(model_to_plot, newdata = status_model_data, re.form = NA)

(status_plot_pdss = status_model_data %>%
  filter(!is.na(target)) %>%
  ggplot(aes(pdss, expected_mean, color = domain, linetype = target)) +
    geom_line(aes(pdss, task_percent, group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_point(aes(pdss, task_percent, group = interaction(subjectID, target, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2, se = FALSE, size = 1.25) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_linetype_manual(name = "", values = c("dotted", "solid")) +
    coord_cartesian(ylim = c(0, 100)) +
    guides(linetype = guide_legend(override.aes = list(color = "black"))) +
    labs(x = "\npubertal stage (PDS score)", y = "social status or academic proficiency\n") +
    dcbw +
    theme(legend.position = c(.8, .14),
          legend.spacing.x = unit(.1, 'cm'),
          legend.box = "horizontal"))

combined

cowplot::plot_grid(status_plot_age, status_plot_pdss, labels = c('A', 'B'), ncol = 2)

SPPC competence analysis

visualize raw

age

(comp_plot_raw_age = comp_model_data %>%
  ggplot(aes(age, competence, color = domain, fill = domain)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_fill_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\npubertal stage (PSD score)", y = "mean competence rating\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

puberty

(comp_plot_raw = comp_model_data %>%
  ggplot(aes(pdss, competence, color = domain)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\npubertal stage (PDS score)", y = "mean competence rating\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

run models

ICC

model.icc = anova(lm(competence ~ subjectID, data = comp_model_data))
round(model.icc$`Sum Sq`[1] / sum(model.icc$`Sum Sq`), 2)

## [1] 0.31

age

`comp model 1` = lmer(competence ~ 1 + domain + 
                              (1 | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 2 age` = lmer(competence ~ 1 + domain + age_c + 
                             (1 + age_c | subjectID), 
                             data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 3 age` = lmer(competence ~ 1 + domain + age_c + age_c2 + 
                             (1 + age_c | subjectID), 
                             data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 4 age` = lmer(competence ~ 1 + domain*age_c + 
                              (1 + age_c | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 5 age` = lmer(competence ~ 1 + domain*age_c + domain*age_c2 + 
                              (1 + age_c | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))

puberty

`comp model 2 puberty` = lmer(competence ~ 1 + domain + pdss_c + 
                             (1 + pdss_c | subjectID), 
                             data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 3 puberty` = lmer(competence ~ 1 + domain + pdss_c + pdss_c2 + 
                             (1 + pdss_c | subjectID), 
                             data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 4 puberty` = lmer(competence ~ 1 + domain*pdss_c + 
                              (1 + pdss_c | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))
`comp model 5 puberty` = lmer(competence ~ 1 + domain*pdss_c + domain*pdss_c2 + 
                              (1 + pdss_c | subjectID), 
                              data = comp_model_data, control = lmerControl(optimizer = "bobyqa"))

compare models

age

anova(`comp model 1`, `comp model 2 age`, `comp model 3 age`, `comp model 4 age`, `comp model 5 age`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = npar,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")

model	model df	AIC	X2	X2 df	p
comp model 1	4	430.67	–	–	–
comp model 2 age	7	429.81	6.85	3	.077
comp model 3 age	8	428.13	3.68	1	.055
comp model 4 age	8	431.13	0	0	–
comp model 5 age	10	431.05	4.08	2	.13

puberty

anova(`comp model 1`, `comp model 2 puberty`, `comp model 3 puberty`, `comp model 4 puberty`, `comp model 5 puberty`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = npar,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")

model	model df	AIC	X2	X2 df	p
comp model 1	4	430.67	–	–	–
comp model 2 puberty	7	431.15	5.51	3	.138
comp model 3 puberty	8	432.13	1.03	1	.311
comp model 4 puberty	8	432.73	0	0	–
comp model 5 puberty	10	435.07	1.65	2	.437

summarize best fitting model

age

# summarize best fitting model
`comp model 3 age` %>%
  broom.mixed::tidy(effects = c("ran_pars", "fixed"), conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  select(-group) %>%
  rename("b" = estimate,
         "SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p)),
         term = gsub("\\(Intercept\\)", "intercept (academic, age 13)", term),
         term = gsub("domainsocial", "domain (social)", term),
         term = gsub("age_c", "age", term),
         term = gsub(":", " x ", term),
         term = gsub("sd__", "", term),
         term = gsub("Observation", "observation", term),
         effect = gsub("ran_pars", "random", effect),
         `b [95% CI]` = ifelse(effect == "fixed", sprintf("%.2f [%.2f, %.2f]", b, conf.low, conf.high), "--"),
         SE = sprintf("%.3f", SE)) %>%
  mutate_if(is.numeric, round, 2) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  select(term, `b [95% CI]`, SE, t, df, p) %>%
  kable(format = "pandoc")

term	b [95% CI]	SE	t	df	p
intercept (academic, age 13)	3.25 [3.16, 3.34]	0.046	71.20	249.20	< .001
domain (social)	-0.22 [-0.30, -0.14]	0.041	-5.38	228.29	< .001
age	0.02 [-0.00, 0.03]	0.010	1.63	73.36	.106
age2	-0.01 [-0.02, 0.00]	0.004	-1.92	315.15	.056

# summarize random effects
print(VarCorr(`comp model 3 age`), comp = c("Variance","Std.Dev."), digits = 4)

##  Groups    Name        Variance Std.Dev. Corr 
##  subjectID (Intercept) 0.022106 0.14868       
##            age_c       0.001837 0.04286  -0.21
##  Residual              0.152041 0.38992

visualize fitted

age

model_to_plot = `comp model 3 age`
REFormString = as.character(findbars(model_to_plot@call$formula)[[1]])
reForm = as.formula(paste0('~(', REFormString[[2]], REFormString[[1]], REFormString[[3]], ')'))
comp_model_data$expected = predict(model_to_plot, newdata = comp_model_data, re.form = reForm)
comp_model_data$expected_mean = predict(model_to_plot, newdata = comp_model_data, re.form = NA)

(comp_plot_age = comp_model_data %>%
  ggplot(aes(age, expected_mean, color = domain)) +
    geom_line(aes(age, competence, group = interaction(subjectID, domain)), alpha = .1) +
    geom_point(aes(age, competence, group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\npubertal stage (PSD score)", y = "SPPC competence rating\n") +
    dcbw +
    theme(legend.position = "none"))

SPPC importance analysis

visualize raw

age

(imp_plot_raw_age = imp_model_data %>%
  filter(!is.na(importance)) %>%
  ggplot(aes(age, importance, color = domain, fill = domain)) +
    geom_point(alpha = .1) +
    geom_line(aes(group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_fill_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    scale_x_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
    labs(x = "\npubertal stage (PSD score)", y = "mean importance rating\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

puberty

(imp_plot_raw = imp_model_data %>%
  filter(!is.na(importance)) %>%
  ggplot(aes(pdss, importance, color = domain)) +
    geom_jitter(alpha = .1, width = .05, height = .05) +
    geom_line(aes(group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", alpha = .2) +
    scale_color_manual(name = "", values = pal_social_academic) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\npubertal stage (PDS score)", y = "mean importance rating\n") +
    dcbw +
    theme(legend.position = c(.88, .2),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

run models

ICC

model.icc = anova(lm(importance ~ subjectID, data = imp_model_data))
round(model.icc$`Sum Sq`[1] / sum(model.icc$`Sum Sq`), 2)

## [1] 0.24

age

`imp model 1` = lmer(importance ~ 1 + domain +
                           (1 + domain | subjectID), 
                           data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))
`imp model 2 age` = lmer(importance ~ 1 + domain + age_c + 
                           (1 + domain | subjectID), 
                           data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))
`imp model 3 age` = lmer(importance ~ 1 + domain*age_c + 
                           (1 + domain | subjectID), 
                           data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))

puberty

`imp model 2 puberty` = lmer(importance ~ 1 + domain + pdss_c + 
                              (1 + domain | subjectID), 
                              data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))
`imp model 3 puberty` = lmer(importance ~ 1 + domain*pdss_c + 
                              (1 + domain | subjectID), 
                              data = imp_model_data, control = lmerControl(optimizer = "bobyqa"))

compare models

age

anova(`imp model 1`, `imp model 2 age`, `imp model 3 age`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = Df,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")

model	npar	AIC	X2	X2 df	p
imp model 1	6	356.28	–	–	–
imp model 2 age	7	356.65	1.62	1	.203
imp model 3 age	8	356.69	1.97	1	.161

puberty

anova(`imp model 1`, `imp model 2 puberty`, `imp model 3 puberty`) %>%
  as.data.frame() %>%
  select(npar, AIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = Df,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")

model	npar	AIC	X2	X2 df	p
imp model 1	6	356.28	–	–	–
imp model 2 puberty	7	357.77	0.51	1	.476
imp model 3 puberty	8	359.27	0.5	1	.482

summarize best fiting model

`imp model 1` %>%
  broom.mixed::tidy(effects = c("ran_pars", "fixed"), conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  select(-group) %>%
  rename("b" = estimate,
         "SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p)),
         term = gsub("\\(Intercept\\)", "intercept (academic, age 13)", term),
         term = gsub("domainsocial", "domain (social)", term),
         term = gsub("age_c", "age", term),
         term = gsub(":", " x ", term),
         term = gsub("sd__", "", term),
         term = gsub("Observation", "observation", term),
         effect = gsub("ran_pars", "random", effect),
         `b [95% CI]` = ifelse(effect == "fixed", sprintf("%.2f [%.2f, %.2f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 2) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  select(term, `b [95% CI]`, SE, t, df, p) %>%
  kable(format = "pandoc")

term	b [95% CI]	SE	t	df	p
intercept (academic, age 13)	3.35 [3.23, 3.48]	0.06	52.25	55.06	< .001
domain (social)	-1.01 [-1.21, -0.81]	0.10	-10.00	55.63	< .001

# summarize random effects
print(VarCorr(`imp model 1`), comp = c("Variance","Std.Dev."), digits = 4)

##  Groups    Name         Variance Std.Dev. Corr 
##  subjectID (Intercept)  0.04644  0.2155        
##            domainsocial 0.19984  0.4470   -0.24
##  Residual               0.29299  0.5413

visualize fitted

imp_model_data_na = imp_model_data %>%
  filter(!is.na(importance))

model_to_plot = `imp model 1`
REFormString = as.character(findbars(model_to_plot@call$formula)[[1]])
reForm = as.formula(paste0('~(', REFormString[[2]], REFormString[[1]], REFormString[[3]], ')'))
imp_model_data_na$expected_mean = predict(model_to_plot, newdata = filter(imp_model_data_na, !is.na(importance)), re.form = NA)

(imp_plot_age = imp_model_data_na %>%
  ggplot(aes(age, expected_mean, color = domain)) +
    geom_line(aes(age, importance, group = interaction(subjectID, domain)), alpha = .1) +
    geom_point(aes(age, importance, group = interaction(subjectID, domain)), alpha = .1) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = .2, se = FALSE) +
    scale_color_manual(name = "", values = pal_social_academic) +
    scale_x_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
    coord_cartesian(ylim = c(1, 4)) +
    labs(x = "\npubertal stage (PSD score)", y = "SPPC importance rating\n") +
    dcbw +
    theme(legend.position = c(.88, .14),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm")))

SPPC competence importance plot

cowplot::plot_grid(comp_plot_age, imp_plot_age, ncol = 2, labels = c("A", "B"))